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^ : 1 1 Abstract 

'^- 2 Analytic and semi- analytic solution are often used by researchers and practicioners to esti- 

'^ ! 3 mate aquifer parameters from unconfined aquifer pumping tests. The non-linearities asso- 

^ [ 4 ciated with unconfined (i.e., water table) aquifer tests makes their analysis more complex 

^ I 5 than confined tests. Although analytical solutions for unconfined flow began in the mid- 

^ ■ 6 1800s with Dupuit, Thiem was possibly the first to use them to estimate aquifer parameters 

'rT) ' '' fro^^ pumping tests in the early 1900s. In the 1950s, Boulton developed the first transient 

>^, 8 well test solution speciahzed to unconfined flow. By the 1970s Neuman had developed solu- 



Oh 



9 tions considering both primary transient storage mechanisms (confined storage and delayed 

10 yield) without non-physical fitting parameters. In the last decade, research into developing 

,__! ' 11 unconfined aquifer test solutions has mostly focused on explicitly coupling the aquifer with 

K*" ■ 12 the linearized vadose zone. Despite the many advanced solution methods available, there 

QQ 13 still exists a need for realism to accurately simulate real-world aquifer tests. 

^' 14 2 Introduction 

o 

^^ 15 Pumping tests are widely used to obtain estimate s of hydrauhc parameters chara c teriz- 

sj ■ 16 ingflow and transport processes in subsurface (e.g., iKruseman and de Ridderl 1990|, iBatu 



% 



17 [1998[). Hydraulic parameter estimates are often used in planning or engineering applica- 

iB tions to predict flow and design of aquifer extraction or recharge systems. During a typical 

19 pumping test in a horizontally extensive aquifer, a well is pumped at constant volumetric 

20 flow rate and head observations are made through time at one or more locations. Pumping 

21 test data are presented as time-drawdown or distance-drawdown curves, which are fltted to 

22 idealized models to estimate aquifer hydraulic properties. For unconfined aquifers properties 

23 of interest include hydraulic conductivity, specific storage, specific yield, and possibly unsat- 

24 urated fiow parameters. When estimating aquifer properties using pumping test drawdown 

25 data, one can use a variety of analytical solutions involving different conceptualizations and 



simplifiying assumptions. Analytical solutions are impacted by their simplifiying assump- 
tions, which limit their applicability to characterize certain types of unconfined aquifers. This 
review presents the historical evolution of the scientific and engineering thoughts concerning 
groundwater flow towards a pumping well in unconfined aquifers (also referred to variously 
as gravity, phreatic, or water table aquifers) from the steady-state solutions of Dupuit to 
the coupled transient saturated-unsaturated solutions. Although it is sometimes necessary 
to use gridded numerical models in highly irregular or heterogeneous systems, here we limit 
our consideration to analytically derived solutions. 

3 Early Well Test Solutions 

3.1 Dupuit 's Steady-State Finite-Domain Solutions 



DupuitI |l857l | considered steady-state radial flow to a well pumping at constant volumetric 



fiowrate Q [ L^/T] in a ho rizontal homogeneous confined aquifer of thickness h [L] . He used 



Darcy's law [Darcyi . Il856| to express the velocity of groundwater flow u [L/T] in terms of 



39 radial hydraulic head gradient {dh/dr) as 

40 where K = kg/u is hydraulic conductivity [L/T], k is formation permeability [L^], g is the 

41 gravitational constant [L/T^], u is fluid kinematic viscosity [L^/T], h = ip + z is hydraulic 

42 head [L], ip is gage pressure head [L], and z is elevation above an arbitrary datum [L]. 

43 Darcy derived a form equivalent to ([T]) for one- dimensional flow through sand-packed pipes. 

44 Dupuit was the first to apply ([1]) to converging flow by combining it with mass conservation 

45 Q = {27rrb) u across a cylindrical shell concentric with the well, leading to 

Q = K{27rrb)^. (2) 

46 Integrating ([2]) between two radial distances ri and r2 from the pumping well, Dupuit eval- 

47 uated the confined steady-state head difference between the two points as 

4B This is the solution for flow to a well at the center of a circular island, where a constant 

49 he ad condition is applied at the edge of the island (r2). 



DupuitI |1857l | also derived a radial flow solution for unconfined aquifers by neglecting the 



vertical flow component. Following a similar approach to confined aquifers, iDupuitI |1857 
estimated the steady-state head difference between two distances from the pumping well for 
unconfined aquifers as 

h'{r,)-h'{n) = ^\og(^]. (4) 
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These two solutions are only strictly valid for finite domains; when applied to domains 
without a physical boundary at r2, the outer radius essentially becomes a fitting parameter. 
The solutions are also used in radially infinite systems under pseudo-static conditions, when 
the shape of the water table does not change with time. 

Equations ([3]) and (|1D are e quivalent when 6 in ([3]) is average head {h{ri) + h{r-2)) /2. 



In developing (jl]), iDupuitI 1857l | used the following assumptions (now commonly called the 



Dupuit assumptions) in context of unconfined aquifers: 

• the aquifer bottom is a horizontal plane; 

• groundwater flow toward the pumping wells is horizontal with no vertical hydraulic 
gradient component; 

• the horizontal component of the hydraulic gradient is constant with depth and equal 
to the water table slope; and 



66 • there is no seepage face at the borehole. 

67 These assumptions are one of the main approaches to simplifying the unconfined flow problem 
6B and making it analytically tractable. In the unconfined flow problem both the head and the 

69 location of the water table are unknowns; the Dupuit assumptions eliminate one of the 

70 unknowns. 
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3.2 Historical Developments after Dupuit 

NarasimhanI |l998) and Ide VriesI |2nn'i 



hydrology and soi l mechanics 



give detailed historical accounts of groundwater 
only history relevant to well test analysis is given here. 
Forchheimerl 1886| first recognized the Laplace equation V^/i = governed two-dimensional 
steady confined groundwater flow (to which ([3]) is a solution), allowing analogies to be drawn 
between groundwater flow and steady-state heat conduction, including the first application 



of conformal mapping to solve a groundwater flow problem. iSlichterl |1898| also arrived at 
the Laplace equation for groundwater flow, and was the flrst to account for a vertical flow 
component. Utilizing Dupuit 's assumptions, iForchheimerl jl898l | develope d the steady-state 
unconflned differential equation (to which (|1]) is a solution), V^/i^ = 0. iBoussinesql |l904 | 
flrst gave the transient version of the conflned groundwater flow equation a^V^/i = dh/dt 
(where a^ = K/Sg is hydraulic diffusivity [L^/T] and Sg is speciflc storage [1/L]), based upon 
analogy with t ransien t heat conduction. 

In Prague, iThienj 1906| was possibly the flrst to use ([3]) for estimating K from pump- 



ing tests with multiple observation wells Simmons! . |2008| . Equation (I3]) (commonly called 
the Thiem equation) was tested in the 1930's both in the fleld (IWenzell |l932| performed 
a 48-hour pumping test with 80 observation wells in Grand Island, Nebraska) and in the 
laboratory ( iWyckoff et al.l 1932| developed a 15-degree unconflned wedge sand tank to sim- 
ulate converging flow). Both found the stea dy-state solution lacking in ability to consistently 



90 estimate aquifer parameters. IWenzell (1942| | developed several complex averaging approaches 



91 (e.g., the "limiting" and "gradient" formulas) to attempt to consistently estimate K using 



92 steady-state confined equations for a finite system from transient unconfined data. iMuskat 



1932l | considered partial-penetration effects in steady-state flow to wells, discussing the na- 
ture of errors associated with assumption of uniform flux a cross the well s creen in a partially 
penetrating well. Muskat's textbook on porous media flow Muskatl . Il937l | summarized much 
of what was known in hydrology and petroleum reservoir engineering around the time of the 
next major advance in well test solutions by Theis. 



3.3 Confined Transient Flow 



TheisI [l935l | utilized the analogy between transient groundwater flow and heat conduction to 
develop an analytical solution for confined transient flow to a pumping well (see Figure [T]) . 
He initially applied his solution to unconfined flow, assuming instantaneous drainage due to 
water table movement. The analytical solution was based on a Green's function heat con- 
duction solution i n an i nfinite axis-symmetric slab due to an instantaneous line heat source 
or sink ICarslaw . 1921 



^ ^ With the aid of mathematician Clarence Lubin, Theis extended 

the heat conduction solution to a continuous source, m otivated to b etter explain the results 
of pumping tests like the 1931 test in Grand Island. iTheid 1935| gave an expression for 



drawdown due to pumping a well at rate Q in a homogeneous, isotropic confined aquifer of 
infinite radial extent as an exponential integral 



(r,t) 
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r^/i^Ost) 



U 



(5) 



109 where s = hQ{r) — h{t, r) is drawdown, Hq is pre-test hydraulic head, T = Kb is transmissivity, 
no and S = Ssb is storativity. Equation ([5]) is a solution to the diffusion equation, with zero- 
in drawdown inital and far-field conditions, 



s{r, t = 0) = s(r -> cx), t) = 0. 



(6) 



112 The pumping well was approximated by a line sink (zero radius), and the source term 

113 assigned there was based upon ([2]), 



, ds 
limr—- 

r-5>o or 
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(7) 



Although the transient governing equation was known through analogy with heat con- 
duction, the transient storage mechanism (analogous to specific heat capacity) was not com- 
pletely understood. Unconfined aquifer tests were known to experience slower drawdown 
than confined tests, due to water supplied by dewatering the zone near th e water table , 
whi ch is re l ated to the formation specific yield (porosity less residual water). iMuskatI 1934 



and iHurstI 1934J derived solutions to confined transient radial flow p roblem s for f inite do- 

I 1 I \i 

mains, but attributed transient storage solely to fluid compressibility. IJacobI |1940l | derived 
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Figure 1: Unconfined well test diagram 



the diffusion equation for groundwater flow in compressible elastic confined aquifers, us- 
ing mass conse rvation and Darcy's law, without recourse to analogy with heat conduction. 



Terzaghil |1923| developed a one- dimensional consolidation theory which only considered the 



125 compressibil it y of the s oil (in his case a clay), unknown at the time to most hydrologists 



Batul . Il998j . iMeinzerl |1928| studied regional drawdown in North Dakota, proposing the 



127 moder n storage me chanism related to both aquifer compaction and the compressiblity of 

12B water. IJacobI |1940| formally showed Sg = Pwdif^p + nf3y^), where py^ and /3w are fluid density 

129 [M/L^] and compressibility [LT^/M], n is dimensionless porosity, and /3p is formation bulk 

130 compressibility. The axis-symmetric diffusion equation in radial coordinates is 



d'^s 1 ds 



a. dt 



131 When deriving analytical expressions, the governing equation is commonly made dimen- 

132 sionless to simplify presentation of results. For flow to a pumping well, it is convenient to use 

133 Lc = 6 as a characteristic length, Tc = Sb'^/T as a characteristic time, and He = Q/(47rT) 

134 as a characteristic head. The dimensionless diffusion equation is 



d'^SD 1 dsD ds 



drjj 



+ 



D 



td Otd dtn 



(9) 



135 where r p = r/Lg, S n = s/Hc, and tu = t/Tc are scaled by characteristic quantities. 



136 The iTheisI |1935l | solution was developed for field application to estimate aquifer hydraulic 



137 properties, but it saw limited use because it was difficult to compute the exponential integral 
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for arbitrary inpu ts. IWerizell |l942| proposed a type-c urve method that enabled graphical 
application of the iTheid J1935J | solution to field data. ICooper and Jacob 1946l | suggested 
for large values of Id {t^, > 25), the infinite integral in the iTheid 1935| solution can be 
approximated as 

f°° e~" fATt\ 

SoitD, rD)= du^ logg ^^ - 7 (10) 

where 7 ~ 0.57722 is the Euler-Mascheroni constant. This leads to Jacob and Cooper's 
straight-line simplification 

Q 



As ^ 2.3- 
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where As is the drawdown over one log-cycle (base 10) of time. The intercept of the straight- 
line approximation is related to S through fflOll This appro ximation made estimating hy- 
draulic parameters much simpler at large to. iHantushl |l961| later extended Theis' confined 
solution for partially penetrating wells. 



3.4 Observed Time- drawdown Curve 



149 Before the time-dependent solution of iTheisI |1935| . distance drawdown was the diagnos- 

150 tic plot for aquifer test data. Detailed dis tance-drawdow n plots require many observation 

151 locations (e.g., the 80 observation well s of We nzel [19365). Re-analyzi ng results o f the un- 

152 confined pumping test in Grand Island, IWenzel [1942.] noticed that the iTheisI 1935| solution 

153 gave inconsistent estimates of Sg and K, attribute d to the delay in the yield of water from 

154 storage as the water table fell. The iTheid |l935j solution corresponds to the Dupuit as- 

155 sumptions for unconfined flow, and can only re-create the a portion of observed unconfined 

156 time-drawdown profiles (either late or early). The effect of the water table must be taken 

157 into account through a boundary condition or source term in the governing equation to 
15B rep roduce observ ed behavior in unconfined pumping tests. 

159 IWaltonI |l960l | recognized three distinct segments characterizing different release mech- 

160 anisms on time-drawdown curve under water table conditions (see Figure |2]). A log-log 

161 time-drawdown plot in an unconfined aquifer has a characteristic shape consisting of a steep 

162 early-time segment, a flatter intermediate segment and a steeper late-time segment. The 

163 early segment behaves like the iTheisI |l935| solution with S = S.^b (water release due to bulk 

164 medium relaxatio n) , the late segment behaves like the lTheisI 1935J solution with S = Sgb+Sy 

165 (Gambolatil . Il976l | (water release due to water table drop), and the intermediate segment rep- 

166 resents a transition between the two. Distance-drawdown plots from unconfined aquifer tests 

167 do not show a similar inflection or change in slope, and do not produce good estimates of 

168 storage parameters. 
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Figure 2: Drawdown data from Cape Co d [Moench et al.l . l200lj . observation well F377-037. 
Upper dashed curve is confined model of lHantushI |l961l | wit h S = S.^b, lowe r dotted curve is 



same with S = Ssb + Sy. Solid curve is unconfined model of lNeumanI |l974| using Sy = 0.23 
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4 Early Unconfined Well Test Solutions 

4.1 Moving Water Table Solutions Without Confined Storage 



The iTheisI |1935j solution for confined aquifers can only reprod uce either the ea rly or late 
segments of the unconfined time-draw down curve ( see Figure [2]). iBoultonI |l954aj suggested 
it is theoretically unsound to use the iTheid jl935| solution for unconfined flow because it 
does not account for vertical flow to the pumping well. He proposed a new mechanism for 
flow towards a fully penetrating pumping well under unconfined conditions. His formulation 
assumed flow is governed by V^s = 0, with transient effects incorporated through the water 
table boundary condition. He treated the water table (where ■?/' = 0, located at z = ^ 
above the base of the aquifer) as a moving material boundary subject to the condition 
h{r, z = ^, t) = z. He considered the water table without recharge to be comprised of a 
constant set of particles, leading to the kinematic boundary condition 



— ih-z) 







(12) 



which is a statement of conservation of mass, for an incompressible fluid. iBoultonI [1954a 
considered the Darcy velocity of the water table as Uz = 



4^^ and Ur 



■^1^, and 



183 expressed the total derivative as 

D _ d 
Di~ di 



Kr dh d K,dh d 



Sy dr dr 



Sy dz dz ' 



184 where K^ and K^ are radial and vertical hydraulic conductivity components. 

185 the kinematic boundary condition f lT2|) in terms of drawdown is 



(13) 
Using (m, 
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at 



Oil 
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Sy dz 



(14) 



BoultonI 1954al | utilized the wellbore and far- field boundary conditions of lTheid |l935| . He 
also considered the aqui fer rests on a n impermeable fiat horizontal boundary dh/dz\^^Q = 0; 



IBS this was also inferred by iTheisI |1935| because of his two-dimensional radial fiow assumption. 



Pagan! |l967l ] extended Boulton's water table solution to the partially penetrating case by 



replacing the wellbore boundary condition with 



, ds 
limrT— 

r^o dr 



b-i<z<b-d 



2nK{i-d) 

otherwise 



(15) 



where £ and d are the upper and lower boundaries of the pumping well screen, as measured 
from the initial top of the aquifer. 

The two sources of non-linearity in the unconfined problem are: 1) the boundary condition 
is applied at the water table, the location of which is unknown a priori; 2) the boundary 
condition applied on the water table includes s^ terms. 

In order to solve this non-linear problem both Boulton and Dagan linearized it by dis- 
regarding second order components in the free-surface boundary condition ( IT^ and forcing 
the free surface to stay at its initial position, yielding 



ds 
di 



K^ds^ 

Sy dz 



ho, 



(16) 



199 which now has no horizontal fiux component after neglecting second-order terms. Equation 

200 flT6l) can be written in non-dimensional form as 



ds 



D 



dto 



-K 



.ds 
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dz 



zd 



(17) 



D 



where K^ = K^/ K^ is the dimensionless anisotropy ratio and a* = S/Sy is the dimensionless 
storage r atio. 

Both iBoultonI |l954a| and iDaganI 1967l | solutions reproduce the intermediate and late 
segments of the typical unconfined time-drawdown curve, but neither of them reproduces 
the early segment of the curve. Both are solutions to the Laplace equation, and therefore 
disregard confined aquifer storage, causing pressure pulses to propagate instantaneously 
through the saturated zone. Both solutions exhibit an instantaneous step-like increase in 
drawdown when pumping starts. 
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4.2 Delayed Yield Unconfined Response 

BoultonI |l954b| extended Theis' transient confined theory to include the effect of delayed 



yield due to m ovement of the water table in unconfined aquifers. Boulton's proposed solutions 
(Il954bl . Il963l ) reproduce all three segments of the unconfined time-drawdown curve. In his 
formulation of delayed yield, he assumed as the water table falls water is released from storage 
(t hrough d r ainage) gra dually, rather than instantaneously as in the free-surface solutions 
of iBoultonI 1954a| and iDaganI 1967l | . This approach yielded an integro-differential flow 
equation in terms of vertically averaged drawdown s* as 

Sds*' 






Ids* 
r dr 



T dt 



aS,i 



OT 



which Boulton linearized by treating T as constant. The term in square brackets is instan- 
taneous confined storage, the term in braces is a convolution integr al representing storage 
released gradually since pumping began, due to water table decline. iBoultonI |l963l | showed 
the time when delayed yiel d effects becora e negligible is approximately equal to -, leading 



to the term "delay index". iPrickettI 1965 1 used this concept and through analysis of large 



amount of field drawdown data with IBoultonI |1963| solution, he established an empirical 
relationship between the delay index and physical aquifer properties. Prickett proposed a 
methodology for estimation o f S, Sy, K, and a of unconfined aquifers by analyzing pumping 



tests with the Boulton 1963 solution. 



Although Boulton's model was able to reproduce all three segment of the unconfined 
time-drawdown curve, it failed to explain the physical mechanis m of the delayed yield process 
be cause of the non-p hysical nature of the "delay index" in the IBoultonI 1963J solution. 



Streltsoval |1972aj developed an approximate solution for t he decline of the water table 



and s* in fully penetrating pumping and observation wells. Like IBoultonI 1954bj , she treated 
the water table as a sharp material boundary, writing the two-dimensional depth-averaged 
flow equation as 

Ids* S fds* d£ 



dh* 



+ 



T \dt 



(19) 



The rate of water table decline was assumed to be linearly proportional to the difference 
between the water table elevation ^ and the vertically averaged head b — s*, 

dt Sybz 

where hz = b/3 is an effective aquifer thickness over which water table recharge is distributed 
into the deep aquifer. Equation (l20l) can be vie wed as an approx ima tion to t he ze ro-order lin- 



{s*-b + 



(20) 



earized free-surface boundary condition (fT6ll of lBoultonI |1954al | and lDaganI |1967l |. Streltsova 
considered the initial condition ^{r,t = 0) = b and used th e same boun dary condition at th e 
pumping well and the outer bo undary (r — >■ oo) used by iTheisI 19351 and IBoultonI |l963| . 



Equation (1191) has the solution [Streltsoval . 11972b 



dt 



-—- = —aT 



d^* 
dr 



(21) 
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where ar = Kz/{Syhz). Substituting (12T1) into (!20l) produces solution (TT8l) of iBoulton 



1954bl . [l963|; the two solutions are equivalent. Boulton's delayed yield theory (like that 
of Streltsova) does not account for flow in unsaturated zone but instead treats water table as 



material boundary moving yertical l y dow nward under influence of gravity. iStreltsoval |1973 



used field data collected by iMeyerl [19621 ] to demonstrate unsaturated flow had virtually no 
impact on the observed delayed process. Although Streltsova's solution related Boulton's 
delay index to phys i cal aq uifer properties, it was later found to be a function of r |Neumanl . 



19751 . iHerrera et al.l . Il978l |. The delayed yield solutions of Boulton and Streltsova do not ac- 
count for vertical flow within the unconfined aquifer through simplifying assumptions; they 
cannot be extended to account for partially penetrating pumping an d obsery a tion w ells. 



Prickett's pumping test in the vicinity of Lawrenceville, Illinois |Prickettl . Il965| showed 



that specific storage in unconfined aquifers can be much greater than typically observed 
values in confined aquifers - possibly due to entrapped air bubbles or poorly consolidated 
shallow sediments. It is clear the elastic properties of unconfined aquifers are too important 
to be disregarded. 
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4.3 Delayed Water Table Unconfined Response 

Boulton's (Il954bl . Il963l ) models encountered conceptual difficult y explaining th e physical 
mechanism of water release from storage in unconfined aquifers. iNeumanI 1972l | presented 
a ph ysically based mathem atic al model that tre ated the unconfined aquifer as compressible 
(hke lBoultonI Il954bl. 1196311 and Streltsoval 1972a]Jbl| ') and the water table as a moving material 



boundary (like iBoultonI 1954a 



and iDaganI 1967|). In Neuman's approach aquifer delayed 
response was caused by physical water table movement, he therefore proposed to replace the 
ph rase "delayed yi eld" by "delayed water table respo nse" . 



Neumaru [1972| | replaced the Laplace equation of iBoultoru |1954al | and iDaganI |1967l | by 



the diffusion equation; in dimensionless form it is 






, 1 dsD , r^ d^SD 



ds 



D 



dt 



D 



(22) 



Like iBoultoru 1954a| and lDaganI 19671], Neuman treated the water table as a moving material 
boundary, linearized i t (using flT7l)\. and treated the anisotropic aquifer as three-dimensional 



axis-symmetric. Like iDaganI [19671 ]. iNeumanI [19741 ] accounted for partial penetration. By 



including confined storage in the governing equation fl22l) . Neuman was able to reproduce 
all three parts of the observed unconfined time-drawdown curve and produce parameter 
estimates (including the ability to estimate Kz) very similar to the delayed yield models. 

Compared to the delay index mode ls, Neurnan's solution produced similar fits to data 
(often underestimating Sy, though), but INeumanI jl975Ul979J question ed the physical n ature 
of Boulton's d elay index. He performed a regression fit between the iBoultonI 1954b| and 
Neumaru 1972| solutions, resulting in the relationship 



a 



Syb 



3.063 -0.567 log 



( Kpr 
V &2 



2\ 1 



(23) 



10 



demonstrating a decreases linearly with logr and is therefore not a characteristic aquifer 
constant, ^ yhen ignor i ng the logarithmic term in (123 p the relationship a = 3Kz/{Syb) 
proposed by IStreltsoval |1972aj is approximately recovered. 

Af ter comparative analysis of various methods for determination of specific yield, iNeumanI 



1987l | concluded the water table response to pumping is a much faster phenomenon than 
dr ainage in t he un saturated zone above it. 



Malamal 201 1| recently proposed an alternative linearization of f ll4p . approximately in- 



cluding the effects of the neglected second-order terms, leading to the alternative water table 
boundary condition of 



^ ds ^^ f ds ^d'^s 



hn 



(24) 



where /3 is a linearization coefficient [L]. The parameter /3 provides additional adjustment 
of the shape of the intermediate portion of the time-drawdown curve (beyond adjustments 
possible with Kd and cr* alone), leading to improved estimates of Sy. When /3 = ( 12^ 
simplifies to flT6|) . 
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4.4 Hybrid Water Table Boundary Condition 

The solution of iNeumanI 19721 . Il974| was accepted by many hydrologists "as the pre- 
ferred model ostens ibly because it appe ars to make the fewest si mplifying assumptions" 
Moench et al.l . l200ll |. Despite acceptance. iNwankwor et al.l 1984l | and lMoenchI |1995] pointed 
out that significant difference might exist between measured and model-predicted draw- 
downs, especially at locations near the water tab le, leading to s ignificantly underestimated 



Sy using Neuman's models (e.g., see Figure E]). iMoenchI |l995j attributed the inability of 



Neuman's models to give reasonable estimates of Sy and capture this observed behavior near 
the water table du e to the later's disregard of "gradual drainage". In an attempt to re- 
solve this problem, iMoenchl 1995| replaced the in stantaneous mov ing water table boundary 
condition used by Neuman with one containing a iBoultonI 1954bj delayed yield convolution 
integral. 



r^ M 

OS sr-^ 



0^7776 



m=l 



Sy dz 



(25) 



This hybrid boundary condit ion ( M = 1 in iMoench 



term 



Boulton 1954b. 1963 and Streltsova 



1972a 



1995| ) included the convolution source 



b| used in their d epth-averaged gov 



Moench 1995 



m- 



erning flow equations. In addition to this new boundary condition, 

eluded a finite radius pump i ng w ell with wellbore storage, conceptually similar to how 



Papadopulos and Cooper Jr.l |1967| modified the s olution o f iTheisI |1935| . In all other re 



spects, his definition of the problem was similar to INeumanI |1974 . 

Moench's solution resulted in improved fits to exp erimental data and produced more 
realistic estimat es of specific yi eld Moench et al.l l200l|, including the use of multiple delay 
parameters am Moenchl . l2003|]. iMoench et al.. [200ll | used ( l25l) with M = 3 to estimate 
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hydraulic parameters in the unconfined aquifer at Cape Cod. They showed that M = 3 
en abled a better fit to the observed drawdown data than obtained by M < 3 or the model 



312 of iNeumaru |1974| . Similar to the parameter a in Boulton's model, the physical meaning of 



am is not clear. 



5 Unconfined Solutions Considering Unsaturated Flow 



As an alternative to linearizing the water table condition of lBoultonI |l954al |. the unsaturated 
zone can be explicitly included. The non-linearity of unsaturated flow is substituted for the 
non-linearity of (1141) . By considering the vadose zone, the water table is internal to the 
domain, rather than a boundary condition. The model-data misfit in Figure [2] at "late 
intermediate" time is one of the motivations for considering the mechanisms of delayed yield 
and the effects of the unsaturated zone. 
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5.1 Unsaturated Flow Without Confined Aquifer Storage 



Kroszynski and Daganl |1975l | were the first to account analytical ly for the effec t of the unsat- 
urated zone on aquifer drawdown. They extended the solution of lDaganl |l967l | by accounting 
for unsaturated flow above the water table. They used Richards' equation for axis-symmetric 
unsaturated flow in a vadose zone of thickness L 



-4IKI)--4HI)--(«4 



^ < z <b + L 



(26) 



where a = b + ipa — h is unsaturated zone drawdown [L], ipa is air-entry pressure head [L], 
< klip) < 1 is dimensionless relative hydraulic conductivity, C{iIj) = dO/dip is the moisture 
retention curve [1/L], and 9 is dimensionless volumetric water content (see inset in Figured]). 
The y assurn e d flow in the underlying saturated zone was governed by the Laplace equation 
(like iDaganl |l967|). The saturated and unsaturated flow equations were coupled through 
interface conditions at the water table expressing continuity of hydraulic heads and normal 
groundwater fluxes. 



a 



Vs ■ n = Vcr ■ n 



e 



(27) 



where n is the unit vector perpendicular to the water table. 



To solve the unsaturat e d flow equation (126|1 . iKroszynski and Daganl |1975| linearized (126|) 
by adopting the iGardnerl 1958| exponential model for the relative hydraulic conductivity, 
k{il)) = e'^°-^'^~''^''\ where Ka is the sorptive number [1/L] (related to pore size). They adopted 



the same exponential form for the moisture capacity model, 6{ip) = e 



f^ki^-M'k) 



where ipk 



is the pressure at which kixjj) = 1, /€„ = k^, and ijja = ipk^ leading to the simplified form 
Cjxjj) = SyK ae'^°-'^'^~'^°-\ In the limit as Kfc = Kq — ^ cxd their solution reduces to that of 



Daganl [1967|. The relationsh ip between pressure head and water content is a step function. 



Kroszynski and Daganl |1975l | took unsaturated flow above the water table into account but 
ignored t he effects of cori fined aquifer stora ge, leading to early-time step-change behavior 
similar to iBoultonI 1954a| and iDaganI 1967 . 
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5.2 Increasingly Realistic Saturated-Unsaturated Well Test Mod- 
els 



Mathias and Butlerl |2006| combined the confined aquifer flow equation ( 122|) witli a one- 
dimensional linearized version of ( 126|) for a vadose zone of finite tliickness. Their water table 
was treated as a fixed boundary with known flow conditions, decoupling the unsaturated and 
saturated solutions at the water table. Although they only considered a one- dimensional 
unsaturate d zone, they included the additional flexibility provided by different exponents 



Ka 7^ Kk). iMathias and Butled |2006| did not consider a partially penetrating well, but they 



did note the possibility of accounti ng for it in pr inciple by incorporating their uncoupled 

which considers a partially penetrating 



drainage function in the solution of iMoenchI |1997 
well of finite radius. 



Tartakovsky and Neumad |2007l | similarly combined the confine d aquifer flow equation 



fl22l) . but with t he original axis-symmetric for m of (|26|) considered by iKroszynski and Pagan 



1975l | . Also like IKroszynski and Dagad [1975 



single exp onent Kg 



their unsaturated zon e was characterized by a 
Kk and referenc e pressure head ijja = ijjk- Unlike iKroszvnski and Pagan 



19751 1 and lMathias and Butlerl 2006l | . iTartakovsky and Neumad 2007l | assumed an infinitely 
thick unsaturated zone. 



Tartakovsky and Neumad |2007| demoristrate d flow in the unsaturated zone is not strictly 



vertical. Numerical simulations by lMoencd |2008| showed groundwater rn ove ment in the cap- 



illary fringe is more horizontal than vertical. IMathias and Butlerl [2006l | and iMoencd |2008 
showed that using the same exponents and reference pressure heads for effec tive saturation 
and relative permeability dec reases model flexibility and unde restimates Sy 



Moench 2008 



predicted an extended form of lTartakovsky and Neumad |2007l | with two separate exponents, 
a finite unsaturated zone, and wellbore storage would likely produce more physically realistic 
estimates of S,,. 



Mishra and Neumad |2010| developed a new generalization of the solution of lTartakovsky and Neuman 



2007l | that characterized relative hydraulic con ductivity and water content using Ka ^ Kk 



'4'a 7^ '4'k and a finitely thick unsaturated zone. iMishra and Neumad 2010l | validated their 
solution against numerical simulations of drawdown i n a s ynthetic aquifer with unsatu- 
rated properties given by the model of Ivan Genuchted 1980 1. They also estimated aquifer 



parameters fro r n Cap e Cod drawdown data [Moench et al.l. l200ll|. comp aring estimated 



van Genuchted | 1980| paramet ers with laboratory values Mace et al.l . 11998 



Mishra and Neumad 2011| further extended their 2010l sol ution to include a finite-diameter 



pumping well with storage. iMishra and Neumad |2010l . l201ll | were the first to estimate non- 
exponential model unsaturat ed aquifer properties f rom pumping test data, by curve-fitting 
the exponential mod el to the Ivan Genuchted |l980l | r aodel. Analyzing pumping te st data of 



Moench et all [200 1| (Cape Cod, Massachusetts) and lNwankwor et all |l984l . Il992| (Borden 



Canada), they estimated unsaturated flow parameters similar to laboratory-estimated values 
for the same soils. 
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6 Future Challenges 

The conceptualization of groundwater flow during unconfined pumping tests has been a 
challenging task that has spurred substantial theoretical research in the field hydrogeology 
for decades. Unconfined flow to a well is non-linear in multiple ways, and the application 
of analytical solutions has required utilization of advanced mathematical tools. There are 
still many additional challenges to be addressed related to unconfined aquifer pumping tests, 
including: 

• Hysteretic effects of unsaturated flow. Different exponents and reference pressures are 
needed during drainage and recharge events, complicating simple superposition needed 
to handle multiple pumping wells, variable pumping rates, or analysis of recovery data. 

• Collecting different data types. Validation of existing models and motivating develop- 
ment of more realistic ones depends on more than just sa t urate d zone head data. Other 



data types include va dose zone water conteri t |Meyerl . Il962| . and h ydrogeophysical 



data like microgravity jDamiata and Led . |2006| or streaming potentials jMalama et al. 
20091 



• Moving water table position. All solutions since iBoultonI |1954a| assume the water 
table is fixed horizontal ^(r, t) = ho during the entire test, even close to the pump- 
ing well where large drawdown is often observed. Iterative numerical solutions can 
accommodate this, but this has not been included in an analytical solution. 

• Physically realistic partial penetration. Well test solutions suffer from the complication 
related to the unknown distribution of flux across the well screen. Commonly, the flux 
distribution is simply assumed constant, but it is known that flux will be higher near 
the ends of the screened interval that are not coincident with the aquifer boundaries. 

• Dynamic water table boundary condition. A large increase in complexity comes from 
explicitly including unsaturated flow in unconfined solutions. The kinematic boundary 
condition expresses mass conservation due to water table decline. Including an anal- 
ogous dynamic boundary condition based on a force balance (capillarity vs. gravity) 
may include sufficient effects of unsaturated flow, without the complexity associated 
with the complete unsaturated zone solution. 

• Heterogeneity. In real-world tests heterogeneity is present at multiple scales. Large- 
scale heterogeneity (e.g., faults or rivers) can sometimes be accounted in analytical 
solutions using the method of images, or other types of superpostion. A stochastic 
approach Neuman et al.l . |2004| could alternatively be developed to estimate random 
unconflned aquifer parameter distribution parameters. 

Despite advances in cons i dering physically realistic unconfined fio w, most real-world un - 
confined tests (e.g., IWenzell |l942i . iNwankwor et al.l |l984j . Il992| . or iMoench et aD |200l| ) 
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419 exhibit non-classical behavior that deviates from the early-intermediate-late behavior pre- 

420 dieted by the models summarized here. We must continue to strive to include physically rel- 

421 evant processes and representatively linearize non-linear phenomena, to better understand, 

422 simulate and predict unconfined flow processes. 
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